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Abstract 


Flows in 3D complex geometries mainly curved ducts with square and circular 
cross sections are studied using a finite volume formulation. The finite volume 
formulation uses a collocated grid arrangement with cartesian components of 
\elocity, and momentum interpolation to avoid pressure- velocity decoupling. 
The numerical computations are carried out for 90^ and 180° bends of square 
cross section and a 90° bend of circular cross section. For a range of curvature 
ratios (R/d) and Reynolds numbers, axial {Utf,) and secondary velocities {Ur) 
are plotted and compared with experimental results of established repute. In 
all cases, the computed and experimental results are close, and validate the 
accuracy and robustness of the Navier Stokes solver. 
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Chapter 1 
Introduction 


The prediction of incompressible fluid flow is needed in almost every branch 
of engineering. Examples include chemical processing plants, nuclear reactors, 
heating and ventilation of rooms, cooling of electronic equipments, meteorol- 
ogy, turbomachines, ground water flows, etc. As the field of computational 
fluid dynamics continues to mature, there is a need of exploiting the most 
recent advances in numerical mathematics, computing architectures and hard- 
ware to solve fluid flow problems of increasing complexity, for industry and in 
fundamental research. 


Successfully obtaining physically meaningful solutions from numerical 
flow predictions depends on (a) the correctness of the mathematical approx- 
imations in the modeling of the physical systems, and (b) the accuracy and 
efficiency of the numerical scheme and its ability to handle the complexity of 
the flow geometries under consideration. 


The choice of an appropriate numerical method often depends upon the 
geometry of the flow domain. The most commonly used solution methodologies 
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are finite-difference, finite-volume and finite-element methods. Both finite- 
volume and finite-element methods belong to the class of weighted residual 
methods[l]. However, in implementation, finite-volume discretizations more 
closely resemble finite-difference methods. Finite-difference methods[l], in gen- 
eral, are algorithmically simplest, while finite-eleraents[2] are the most flexible 
with regards to flow geometry, but are rather involved algorithmically. The 
finite- volume method is often seen as a good compromise between the two - 
being simpler to implement than the finite-element method and being more 
suitable for complex domains than the finite-difference method. The finite- 
volume method has the additional advantages of discretizing the conservation 
form of the governing equations. This implies that the descretised equations 
preserve the conservation laws. This is a particular advantage w'hen obtaining 
accurate solutions for internal flows or flows with shocks. 


If the geometry is rectangular then finite-difference methods are often 
preferred for their simplicity. However, with the growing use of CFD in in- 
dustrial applications, computations are done for increasingly complex geome- 
tries. The recent advances in grid generation techniques[3, 4, 5] allow the use 
of boundary-fitted grids which permit selective refinement near the boundary 
walls, etc., and thus reduce storage, computational time, and improve accuracy. 
However finite-difference schemes used on body-fitted coordinates require us to 
map the complex domain onto a regular rectangular grid in which transformed 
equations are solved. This could be global mapping, i.e., mapping the entire 
flow' domain onto a simpler computational domain, or, more commonly, local 
mapping, i.e., mapping each successive grid point and its related neighbours 
onto a succession of small and simple domains. The transformation approach 
is often adopted even for finite-volume schemes. 


The most commonly used grids for finite-difference and finite-volume 
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methods are either staggered or non-staggered (collocated) grids. The stag- 
gered grid arrangement is quiet robust and efficient in orthogonal grid compu- 
tations and avoids pressure-velocity decoupling associated with non-staggered 
grids. On uniform grids, the grid points have a simple relationship with the 
grid index and their locations need not be stored. However, on non-uniform 
grids, the coordinates of each grid point need to be stored. On such grids, a 
staggered system will store many times more information than non-staggered 
ones, as each component of velocity and the pressure have different locations in 
the former. In addition, the grid transformation parameters for each variable 
is different on staggered grids, and again need to be stored. This also poses 
problems in the discretization of the non-linear terms on a non-uniform mesh. 
All this is heavy price to pay for the staggered grid’s pre-eminent advantage 
— the avoidance of pressure-velocity decoupling. 


On the other hand, the collocated grid arrangement stores all the depen- 
dent variables at the same grid points and only one set of grid relations needs 
to be considered. This arrangement is susceptible to pressure-velocity decou- 
pling, but this can be avoided by using momentum interpolation [6]. In view 
of the above considerations, the collocated grid arrangement seems preferable 
for fluid flow calculations in complex geometries. 


In this thesis we attempt to solve incompressible 3-D flows in complex 
geometries using finite volume formulation that was developed for DRDL, Hy- 
derabad, as part of a project. It has been thoroughly tested for 2-D and 3-D 
regular geometries. This thesis validates the code for complex 3-D geome- 
tries. The finite volume formulation is described in the next chapter. The 
scheme uses a pressure based primitive variable formulation and general hex- 
ahedral elements which can be used to represent complicated computational 
domains. A collocated grid arrangement is chosen for the storage of the de- 
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pendent variables to facilitate calculation in complex domains. The problem 
of pressure- velocity decoupling is avoided using momentum interpolation. The 
Cartesian components of velocity and pressure are used as dependent variables. 
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Chapter 2 


Finite - Volume Formulation 


The finite volume formulation presented below has been developed by 
DrA^Eswaran and co-workers as part of a project for DRDL. Hyderabad. The 
full details are given in the project report [7]. 


2.1 Governing Equations 


The three-dimensional Navier-Stokes equations for laminar flow for an arbi- 
trary spatial region of control volume bounded by a closed surface S can be 
expressed in the following general convection-diffusion-source integral form: 

■^ / pd\’ + f pu • dS = 0 (2.1) 

ot Jv Js 

f p(f>d\' + [ [pud - • dS = f S^dV (2.2) 

ot Jv Js Jv 

where p represents the fluid density, u is the fluid velocity, 0 stands for any 

vector component or scalar quantity, Sp is the volumetric source term. 

For incompressible flow of Newtonian fluid, the equations takes the form 

^u-dS = 0 ' (2.3) 
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Description of the Finite-Volume 


+ = (2.4) 

and the source terra for momentum equations becomes - • dS where 

I is the unit tensor. 

In this formulation we work with cartesian components of velocity. 
So 4> can be the three cartesian components of velocity u.v.w, as well as any 
scalar e.g., temperature, which needs to be determined. 


2.2 Description of the Finite- Volume 

Figure 1 shows the kind of control volume used, with eight vertices but other- 
wise irregular. 

The conservation equations are discretized by employing the general 
finite-volume approach. The solution domain is subdivided into a number 
of contiguous (finite) control volumes. The control volumes are defined by the 
coordinates of their vertices which are assumed to be connected by straight 
lines. The coordinates of the control volume vertices are calculated by some 
grid generation procedure. A Collocated Grid arrangement is employed, 
i.e., all the dependent variables {u,i\w,p) are defined at the same location, 
the centroid of the control volume in Figure 1. The six neighbouring control 
volume centers are indicated by E, W, N, S, T, and B (for the east, west, north, 
south, top and bottom neighbours). The face center points e, w, n, s, t, and b 
are located at the corresponding face centroids of control volume P. Note the 
edge-centers te, be. ne, se, etc. which are also needed in the computations. 

Before the time-stepping can proceed all geometrical parameters need to 
be computed as part of the investigation. These parameters include the surface 
vectors for each of the six faces e, w', n, s, t, b of each finite volume and also 
its volume. These computations can be done by the methods shown below. 
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Description of the Finite-Volume 


2.2.1 Surface areas and Volume 

The finite volume vertices are numbered 1 to 8 in the manner shown in figure 
1. The outward surface normals and volume can be found in the following 
manner suggested by Kordulla and Vinokur[8]. 

Defining = r, — Tj where r^, Tj are the position vectors of points i and 
j respectively, we have for 

• East face 


2(^74 X rgs) 


(2.5) 


• West face 


Su- = 2(^16 X rss) 


(2.6) 


• North face 


S„ = -(r27 X res) 


(2.7) 


• South face 


Ss = -(ri8 X r45) 


( 2 . 8 ) 


• Top face 


S,= 


1 


(^75 X res) 


(2.9) 


• Bottom face 
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Discretization Procedure 


S6 = ^(ri3 X rs^) (2.10) 

The volume of the control cell is calculated from the cell corner coordi- 
nates assuming that the cell corners are joined by linear segments to form the 
six cell faces. 

V=^Tn-{Ss + St + S,,) ( 2 . 11 ) 

2.3 Discretization Procedure 

All the conservation equations have the same general form, represented by 

podV + Ijpn4> - • dS = ^ S^dV (2.12) 

The main steps of the discretization procedure to calculate convection and 
diffusion fluxes and source terms are outlined below. The rate of change and 
source terms are integrated over the cell volume, whereas the convection and 
diffusion terms form the sum of fluxes through the CV faces. 

2.3.1 Continuity Equation 

Equation (2.3) is discretized in the followung way. 

J pu dS^ p{u- S)j = puj • Sj (2.13) 

^ j=€,w,n,s,t,b j 

w'here Sj is the surface vector representing the area of the cell face and Uj 
is the velocity defined at the face center j. 

In discretized form the continuity equation follow's 

YFj = Fe + F^. + Fr, + F, + Ft + Ft^(i (2.14) 

J 

w'here the Fj is the outward mass flux through face j, defined by 

Fj = puj • Sj 
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Finite - Volume Formulation 


Discretization Procedure 


2.3.2 Discretization of General Equation 

(a) Rate of change 

The value of the dependent variable <p at the centroid of the control volume 
(the geometric center) represents an average over the C\' as a whole. Thus 

where V is the volume of the cell. 


(b) Convection fluxes 

The surface integral over convection flux of variable (j) can be approximated in 
the following form 


f pu0-dS Si ^pOj(u-S)j 
3 

where Oj is the value of (p at the center of face j. Thus 

pU(;^.dS Fe(l>e + Fw't^w + Fnd„ 4- Fs(l>s + Ft(pt + 


( 2 . 16 ) 


( 2 . 17 ) 


where de is the (interpolated) value of the variable (p at the east face center, 
etc. This can be evaluated by using a ‘center difference’ linear interpolation 
between the neighbouring nodal values Op and (pp- At east face the value of 
(pe is given by 


— 


1 E 


TT^P + 


V 


—(pE 


( 2 . 18 ) 


Vs-i-lp Te + Tp 
where Ve and Vp are volumes of the cells around the points E and P respec- 
tively and (pB and (pp are the values of the dependent variable at these points. 
Since all dependent variables are defined at the same location, exactly the same 
interpolation scheme is used to express all of them at the interfaces. However 
to use the center difference approximation above to compute the convection 
fluxes may lead to numerical stability problems. Therefore the convection flux 
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Discretization Procedure 


is split into a first-order upwind differencing scheme (UDS). and another part, 
which equals the difference between the central difference scheme (CDS) and 
UDS approximations; 

F.Oe = + 7[iFe<Pef^^ - {Fe(pef^^] (2-19) 

The upw'ind differencing scheme is based on the assumption that the convected 
cell face value is equal to that at the upstream cell along the same coordinate 
direction. Thus, the value Oe at the east face is assigned the value cpp if 
Ue > 0, i.e., the flux F^ is positive, and the value (pp if Ug < 0, i.e., the flux Fg 
is negative. This can be conveniently summarized as 

Fe0e = M[FeM-OE[hFg,0]] + j{Fg{-^^4>p + -A^cl>E) 

- <^p[[F„0]]^Oe[[-F„0]]} (2.20) 

Here [[p, 9]] denotes the maximum of p and q. Similar expressions can be 
written for rest of the faces 

Fu-(pw = <f>p[[Fu: 0]] - Oh-[[— F^ y, 0]] -f 7{i'V'(r; — <P vv + rz — ^ 

v'lv + t p ^\v -r V p 

- «j!>p[[Fu.. 0]J -r Ou-[[-Fu,, 0]]} (2.21) 

Fn^n = 0p[[Fn.O]] - O.v[[-F„.0]] -f- 7{F„( _. r <^P T .. 

V > + V A' \ p ^ y N 

- .#f[|F„.0)] + o,vl[-F„,0)I} ^ _ (2,22) 

FA. = »j.[[F.. 01] - Osll-fi. 0)1 7{C(- . 7 

V 5 + 1 P V S + ‘ P 

- «Sp[|C.0ll-K>s|(-C,0]]} _ ^ (2.23) 

FA, = Mfi-Oll-<>rll-S,O]l^-7{F,(T7-Ur^0p-HYTl£Y^0T) 

yp -j- V p \ p ’T vp 

- ?!*p[[Ft.0]] + c>r[[-Ft,0]]} (2.24) 

Fb4>b = <l>p[[Fb, 0]] - ob[[--^6) 0]] + 'y{Fb{Tr~TT^^B + p~~7Tr^p) 

- <;ip[[F6.0]] + (^B[h-F6,0]l} (2.25) 

If a fully implicit method is used for time-stepping, the upwind parts of the 
above equations are implicit' in that they are incorporated in the coefficients of 

“ 'lENTRAL UBRAW 
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Finite - Volume Formulation 


Discretization Procedure 


the unknown velocity during the pressure-velocity iterations. The CDS terms 
on the other hand are evaluated using previous iteration values and used as 
source term on the right hand side of the same equation. This is the so-called 
‘deferred correction’ approach of Khosla and Rubin[9]. Multiplication of the 
explicit part by a factor ') {0 < 7 < 1) allows the introduction of numerical 
diffusion. (7 = 0 means pure UDS, 7 = 1 means pure CDS). The deferred cor- 
rection approach enhances the diagonal dominance of the coefficient matrix, 
which adds to the stability of the solution algorithm. 


(c) Diffusion fluxes 

The diffusion flux of variable 0 through the cell faces can be evaluated as 
follows 

I r^Vd> • dS ^ (r^Vd> ■ s),- = Y. -Ff ( 2 - 26 ) 

^ j=e,u'.n,s,t,6 j 

For any face we can write 

Sj = Qjn' -f Q2n^ -I- Q'an^ (2.27) 

where n^, and are any three linearly independent (not necessarily or- 
thogonal) unit vectors. Therefore 

V^-Sj = W-(Qin'-f Qsn^ + aan^) (2.28) 

= QiVd • -f a2V<6 • n^ -t- oaVp ■ 

If A(f>^ , A<fF , A(f>^ are the differences in (f) between the two ends of the line 
segments Ax^, Ax^, then 

Atf^ = V(?!>.Ax\ AdP^ = V(^.Ax^ AcjP = Vo.Ax^ (2.29) 


If furthermore Ax^ Ax^, Ax^ are in the directions of nT n^ and respectively, 

then it follows from equation (2.29) that 


AtjP 

Ax^ 


= Vcj>.n\ 


Ax^ 


= 


A# 

Ax3 


V.?.n^ 


(2.30) 
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Discretization Procedure 


where etc. are the magnitudes of Ax’ etc. Consequently we can write, 
using (2.29) and (2.30) 


„ , ^ Ao’ A02 ^^3 

V<p.S — Qi v~r + (^2 - — T + 


Ax’ 


■ Ax^ 


Ax^ 


To get Q'l, Q 2 and Q 3 , we express 

n’ = (nil "12 7^13) 
= (n2i n22 ^ 23 ) 
= (7731 7232 7X33) 


(2.31) 


(2.32) 


where rxn, 7 x 12 , 72,13 are the cartesian components of n and which can be easily 
determined by where AxJ, Axj, Axg are the three components 

of vector Ax’, etc. The other values n 2 i ..7233 can be similarly determined. 
Therefore, equation 2.27 can be written as 


’ 72,1 

72,2 

7213 ' 


f ai ] 


j 

7221 

7222 

7223 

< 

^2 

^ = < 

522 

7231 

7232 

7233 . 


1 ^3 J 


1 532 J 


(2.33) 


where Sij, 52 j, Saj are the cartesian components of the surface vector S_,. 


Using Cramer’s rule 


ai 


Di _ D 2 

D' D' 



(2.34) 


where D is the determinant of the coefficient matrix. Di is obtained by re- 
placing the first column of D by the column with elements Sij,S 2 j.S 3 j. etc., 
we can get ai, 02 , Q 3 . 


The diffusion flux is made of two distinct parts: normal diffusion and 
cross-derivative diffusion. The second part arises from the nonorthogonality of 
the grid. The normal derivative diffusion flux of 0 through any cell face involves 
the values of 4> at cell centers whereas the cross-derivative diffusion flux takes 
into account the edge center values of <!>. (The normal derivative diffusion flux 
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Discretiza tion Proced ure 


is treated implicitly and is coupled with the implicit part of the convective 
flux to calculate the main coefficients of the discretized equations while the 
cros-^-derivative diffusion flux is treated explicitly to avoid the possibility of 
producing negative coefficients in an implicit treatment. This term together 
with explicit part of convective flux is added to the source term). The example 
of the east face is taken to illustrate the diffusion model. Given the edge center 
values o-f. Ofeg. d’se> d'ne ''‘6 can get the normal diffusion term ^ , and the 
cross diffusion term and . We can find ai, 02 . 03 by the procedure 

outlined above and finally compute the diffusion flux by 




-r<p(Q 


0E 


Ax' 


Op Of 

+ 02 — 


Ax= 


+ <^3' 




^be 




(2.35) 



Figure 2.2: Face Reitresentation 

To calculate the edge center values appearing in cross-derivative diffusion 
flux the following interpolation scheme is proposed. 


, ^ TE ^ P ^ P 

(f)ie = Tt— ^P ^ —On - —Oe 


1 


-Or 


{2.3Gi 


tot 


' ioi 


' tot 


tot 
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Discretization Procedure 



I tot 


ITe + lp + lr-l-l£; 


II 

I BE 

^ tot 

+ 

TT—QBE + TT— Oe + T7—<Pb 

* iot ^ tot *'tot 

(2.37) 


I tot 

zn 

I BE + If+1b'^I£ 


Qne — 

I'a’e 

r . Qp 

*tot 

+ 

Fp . 1 .V , Ve j 

rr-<pNE + tc-Oe + 7j~d>N 

tot * tot ^tot 

(2.38) 


Vtot 


I A'E + I'P + 1 .V + I'E 



VsE . 

r ^ 

^ tot 

+ 

t^4>se + y^oe + 

^ tot ^ tot ^ tot 

(2.39) 


^ tot 


1 'SE + I "p + ^ i + I 'e 



where Vte is the volume of the cell to the top-east of cell P etc. Other edge 
centers can be similarly interpolated. 

(d) Sources 

The source term is to be integrated over the cell volume. By assuming that 
the specific source at the ex'" center represents the mean value over the whole 
control volume, we can write 

j^S^dV^{S^)pV (2.40) 

Apart from the real source 5^, explicitly treated parts of the convection and 
diffusion fluxes may also be added to The momentum equation contain an 
extra term (pressure). This term is also treated explicitly. Its discretization 
is analogous to that of the ordinary diffusion flux, i.e., for the momentum 
equation the pressure term is 

- / pn^S « -Y^^PjSrj 

where pj is the pressure at the j'* face center and S,j is the component of 
the surface vector for face j. 
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Semi-Explicit Method 


2.4 Semi- Explicit Method 

Existing methods which solve the incompressible Xavier-Stokes equations fall 
into two categories, namely, semi-explicit and implicit schemes. In the first, the 
momentum equations are discretized in an explicit manner with the exception 
of the pressure gradient terms, which are treated implicitly; the continuity 
equation is also enforced implicitly. As a consequence, the equation-coupling 
reduces to a Poisson equation for the pressure corrections. Such schemes, 
however, because of their reliance on explicit differences, suffer from time- 
step restrictions. In the second category, the equations are discretized fully 
implicitly, with the coupling being more complicated. Some methods which 
fall within this group emplo}* sequential iteration (such as the SIMPLE method 
and SIMPLER), in which the equations for each variable are solved repeatedly 
in succession. 

2.4.1 Semi-Explicit Time-Stepping 

For the present situation we will adopt semi-implicit scheme in w’hich the 
equation 

+ Fr = - (2-41) 

has to be solved along wdth 

Y,Fp'=0 (2.42) 

j 

for each finite volume cell. We adopt a two step process. First a predicted 
velocity u* is found which satisfies the equation 

= -EpPd (2-43) 

J J 

Subtracting equation (2.43) from (2.41) we get 
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Pressure - Velocity Corrections 


where the corrections 

(2.45) 

(2.46) 

2.5 Pressure - Velocity Corrections 
(Pressure Correction Poisson Equation) 


u'p = — Up 




Pj=Pj -Pj 

and the corresponding flux corrections F^ have to satisfy 

yF; = -yF; 


Due to the non-staggered variable arrangement, if the variables (velocities and 
pressure) at the cell faces are calculated by linear interpolation between the 
adjacent cell centered quantities then the pressure-velocity iterations do not 
converge and lead to a checker board pressure field. Therefore it is important 
to use momentum interpolation[6] in which the velocity at the cell faces are 
computed by allowing linear interpolation of the convective and diffusive terms 
but not of the pressure term. This is done in the following way: 

If the discretized equation for velocity is written as 

up = uf - ^(F^ -H F^) + (2.47) 

P\ p p\p 

where Up'* is the n''" time-step velocity and is the pressure term. We define 
a new variable vp which is 

vp==uf--^{F^p + F^) (2.48) 

pVp 

Now to estimate the mass flux Fe — pUg • Sg at the east face, straight forward 
linear interpolation would use. 

Fe = p(up, Up) • Sg (2.49) 
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where the over bar indicates a linear interpolation using equation (2.18). But 
momentum interpolation would use 


Fe = jO(vp. Vf) • Se - AtVp • Se (2.50) 

where the gradient Vp is estimated using pp,pjs and also other neighbours for 
a non-orthogonal grid. 

The velocity and pressure fields are calculated with the following Gauss- 
Seidel type algorithm; 

1. Use equation (2.47) to compute the cell-center Up and vp. Make an 
initial guess pp. Use interpolation (Eq. 2.18) to obtain the face-center 
quantities vj and Pj. 

2. Compute the mass flux through each cell face j using 

F; = pVj-Sj-AtVpj-S, (2.51) 

3. Use equation u' = — y Vp' to compute the flux correction at the face j, 

F; = -MVp'j ■ S, (2.52) 

This is computed using the formulation for diffusion fluxes [equations 
(2.28-2.35)] replacing (p by p'. 

4. Compute residual for each cell. 

s = (2'“) 

j 0 

5. Calculate the cell-center pressure correction from the relation 

p'p —p'p + uj— (2.54) 
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where is relaxation factor and op stands for diagonal coefficients that 
is calculated from 


6 . 

7. 


ap = —At 


ai 


<^2 I . Qs 

Axi ^ Aj2 Ax 3 Ax3 


Cll , 


(2.55) 


where ai,Q 2 , etc. are the same as in equations (2.28-2.35). 
If ^rms > e goto step 3 

Store the updated mass flux through cell-faces from 


F* = F* - AfVp; ■ Sj 

8. Store the updated pressure at cell-center pp = pp V p'p 

9. Store the cell-center corrected velocity 

Mp - 

Up = Up -t- Up 


(2.56) 


(2.57) 

(2.58) 


10. Goto step 1 and repeat the process until steady state is reached. 


It can be shown that satisfying equation (2.46) is the same as solving 




Vu* 

At 

(2.59) 

and finding the corrections 


u' = 

-AtVp’ 

(2.60) 

such that 


Un+l 

= u* -f u' 

(2.61) 


pT,+l 

= p"+p' 

(2.62) 


Thus the procedure above can be called the pressure correction Poisson 
solver. 
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Initial and Boundary Conditions 


It should be noted that in this code we compute a pseudo pressure field 
which allows us to satisfy the continuity equation. The pseudo pressure (or 
pressure correction approach) allows us to use simple boundary conditions for 
pressure. 


2.6 Initial and Boundary Conditions 


The following boundary conditions can be incorporated; 


• No-slip boundary conditions 


u = 0 


• Specified boundary conditions 


U = Uo 


• Free-slip boundary conditions 

u„ = 0, 


(9uf 

dn 


= 0 


where n represents the normal to the boundary. The boundary conditions for 
the pressure correction are 


For pressure correction = 0 

The boundary conditions for velocity and pressure correction are implemented 
easily using fictitious cells. 
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Results and Discussion 


3.1 Flow In Curved Ducts 

Flows in curved ducts occur in diverse applications such as centrifugal pumps, 
aircraft intakes, river bends and in the cooling coils of heat exchangers. A 
characteristic which distinguishes such flows from those in straight ducts is 
(a) the generation of stream wise vorticity or "secondary motion” within the 
duct, resulting in a pressure loss, (b) the spatial redistribution of stream wise 
velocity and (c) increased heat transfer at the duct walls. The nature of the 
secondary flow depends, among other parameters, on the curvature ratio of 
the pipe 5=R/d (where R is the mean radius of curvature of the bend and d 
is the hydraulic diameter of the duct) and the flow Reynolds number. 


The present work consists of numerical experiments on 90° and 180° 
bends of square and circular cross section with strong curvature. For the 
problems considered the steady state solution is obtained by the false-transient 
approach, i.e., by time marching the unsteady governing equations. The un- 
steady Navier-Stokes equations describing the flow are : 
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(3.1) 

(3.2) 

(3.3) 

(3.4) 


In the above equations, the velocities have been nondimensionalized with 
respect to the inlet velocity U, and all lengths have been nondimensionalized 
with respect to the hydraulic diameter d. The Reynolds number, is given by 
Re = where Vc is the bulk mean velocitv and i/ is the kinematic viscosity 
of the fluid. 


3.2 Flow In A 90° Bend Square Curved Duct 

3.2.1 Problem Description 


The problem geometry is shown in isometric view in the figure 3.1 . It consists 
of a 90° bend of unit square cross section with a mean radius of curvature of 
R=2.3 and with straight extensions upstream and downstream of lengths that 
are 45 times the hydraulic diameter. 


The flow parameters for this case are Re of 790, based on inflow conditions 
and the hydraulic diameter d (= 4.0XArea/Perimeter) of the duct, curvature 
ratio (5'=2.3(= R/d) where R is the mean radius of curvature of the bend duct. 
The computational domain consists of 64 stations in the stream wdse direction 
(with 30 stations in the bend) and a 24X24 grid in the cross plane. 
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Figure 3.1: Schematic of 90° bend square duct 


The initial conditions, which are arbitrary for this steady state problem, 
are zero velocity conditions. The boundary conditions for the problem are 
given by no slip conditions on the walls and 

u = l,tJ = re = 0 at inlet 

dv 

u = — = It; = 0 at outlet, 

ay 


3.2.2 Results 

In the bend the fluid moves away from the inner radius wall towards the outer 
radius wall under the influence of centrifugal forces. The velocity vector profiles 
in the x-y plane (the plane of symmetry) depicted in the figure 3.2 show' this 
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phenomenon. The prohle is parabolic at the entry of the bend and we can 
observe between 30° and 90° stations a large deceleration in the flow at the 
inner radius and a corresponding acceleration at the outer radius wall. 



Figure 3.2: Velocity vectors in the symmetry plane of a 90° bend duct 


h. 


% • 


1 ; 



The structure of the streamwise velocity is shown in the contours of 
constant azimuthal velocity( U^) in the figures 3.3 and 3.4. The figures demon- 
strate the complexity of the flow. The velocity contours show that with the 
progressive movement, the high velocity region is found nearer to the outer 
wall and along the side walls with a corresponding low velocity region adja- 
cent to the inner w'all. This problem has been solved for the same configuration 
and the same parameter set, as the experiment of Humphrey et al.[10], w’ho 
also numerically simulated their experiment. The figure 3.5 shows the present 
results in comparison w'ith the experimental and computed data of ref[10]. 
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Flow In A 90° Bend Square Curved Duct 


The secondary flow in the bend is shown at different locations by the 
cross flow velocity vectors in the figures 3.6 and 3.7. Figures (a) and (b) show 
the cross sectional flow at (f)=-{F and 90° respectively, while (c) and (d) show 
the same at distances 9 d and 15.75 d downstream of (the 90° point of) the 
bend. It can be seen that the secondary flow is quiet complex and that it 
persists for some distance downstream of the bend. 



Figure 3.8: Measured(+) and computed friction coefficients at the center-line 
of the inside and outside curved walls of the bend. 

The figure 3.8 shows the computed friction coefficients comparing with 
the measured results by Ward-Smith(1971). The coefficients at the outer and 
inner duct walls are plotted as positive and negative on the figure3.8, merelj 
to distinguish between them. While the initial magnitudes are the same, the 
Cf at the outer wall increases dramatically in the bend and does not regain 
fully-developed values till some distance downstream. The Cj at the inner 
w'all does not show such marked variation. The differences arise due to the 
secondary flow structure seen in the previous figures. 
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3.3 Flow In A 180° Bend Square Curved Duct 

3.3.1 Problem Description 

The problem considered is the same in configuration and parameters as the 
experimental work of Phille et al.[14] and is shown in isometric view in figure 
3.9. The model consists of a 180 degree bend of unit square cross section with 
a mean radius of curvature of i2=6.45, and has straight extensions upstream 
and downstream that are 56 times the hydraulic diameter. The inlet section 
is sufficiently longer to ensure fully developed flow. 




Figure 3.9: Schematic of 180° bend square duct 
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The flow parameters are Re of 574 (Re = ^). where V'e is the bulk 
mean velocity and u is the kinematic velocity of the fluid, and a curvature 
ratio R/d)=6.45. This is a very severe test case where the flow is quiet 
complex. The pre<"ent computations are done on a relatively coarse grid with 
60 stations in the bend and a (24X24) grid in the cross plane. 


The initial conditions, which are arbitrary, for this steady state problem, 
are zero velocity conditions. The boundary conditions are no slip conditions 
at the walls and 


n = l,r = t(: = 0 at inlet 

On 

— = V — w = Q at outlet, 

dx 


3.3.2 Results 

Centrifugal forces arising due to the presence of longitudinal curvature causes 
the peak of the axial velocity profiles to shift towards the outer wall. The 
streamwise flow development is summarized in the figure 3.10. 


The flow development in the bend i.s depicted in the figures 3.11 to 3.14 
which show the azimuthal velocity component I'c measured as a function of 
radial coordinate in the mid plane of the cross section. 2=0.5, at intervals of 
18°. The figures show the comparison of the present results with the computed 
results of Kaushik and Rubin ref[13|, obtained on the same sized grid and the 
experimental data[14]. The results are in close agreement with each other. 
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Figure 3.10: Streamwise velocity development in the bend shown as velocity 
vectors plotted in the x-y plane of symmetry. 
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Flow In A 180° Bend Square Curved Duct 


The figure 3.15 shows the presence of two pairs of counter rotating pri- 
mary and secondary vortices. This is in agreement with the experimental data 
which shows the presence of a second vortex pair in the region between 0=108° 
and 171°. 



Figure 3.15: Two dimensional vector profiles of the secondary velocities at 
</)=162°. 
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3.4 Flow In A 90® Bend Circular Duct 

3.4.1 Problem Description 

The flow in a 90° bend of circular cross section is again a stringent test for 
a solver for flows in complex geometry. The problem considered is same as 
in the experimental work of Bovendeerd[15] and is shown in fig.3. 16(a). The 
model consists of a 90° bend in circular cross section pipe of diameter of 8mm 
and a bend of a mean radius of curvature 24mm. The model has a straight 
extension downstream that is 12.5 times the hydraulic diameter. A fully de- 
veloped parabolic profile is introduced at the entrance of the bend (in the 
experimental[15] a long inlet section ensured near parabolic profile at the bend 
entrance). 

The flow parameters are a Reynolds number Re{Re = =700 where 

Vc is the bulk mean velocity and u is the kinematic viscosity of the fluid, and a 
curvature ratio 5{= R/d)=Z. The grid employed for the circular cross section 
is shown in fig.3. 16(b). The grid was generated using a Laplace method[16]. 
[First a circular grid with a diameter of 7.75mm was generated using Laplace 
method which mapped it onto a square. Another grid point was added on each 
side with a diameter of 8.25mm which ensured that the cell points on the outer 
control volumes w'ere exactly at 8mm which represents the physical domain. 
However one cell point at every 45° still could not fall on the ph^^sical domain 
as the control volume is not exactly orthogonal. For this case the boundary 
conditions are imposed on the cell points which represents the physical domain 
and hence no fictitious cells. The zero pressure conditions at the outlet which 
were implemented for the earlier problems created numerical problems and the 
results were obtained by implementing Neumann boundary conditions at the 
outlet]. The computational domain consists of 51 stations in the streamwise 
direction(with 30 stations in the bend) and a 22X22 grid in the cross plane. 
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Flow In A 90° Bend Circular Duct 


The initial conditions for this steady state problem are positive and negative 
parabolic profiles for n and v respectively, and u'=0. The boundary conditions 
are zero velocities at the cells on the cross section, a positive parabolic profile 
of u [with V and w being zero] at the inlet and a Neumann boundary condition 
for V [with u and w being zero] at the outlet. 


3.4.2 Results 

The velocity measurements are shown at seven axial stations in the bend 
corresponding to 4)^ 0°, 4.5°, 12°, 22.5°, 40.5°, 58.5° and 82.5°. [Some of the 
results were interpolated so as to have the results in near comparison to 
the experimental[15] whose axial stations in the bend were at (^=0°,4.6°, 
11.7°, 23.4°, 39.8°, 58.5° and 81.9°]. The development of the axial flow in the 
plane of symmetry is shown in fig.3.17. The curvature has no effect on the 
velocity profiles between 4)—^^ to 12°. Further downstream the maximum of 
the axial velocity([,^,;j,) has shifted towards the outer bend reflecting the strong 
curvature effects. 

0 2 \’ 


Outside of Bend 
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A clear picture of the axial flow is given in the figures 3.18 to 3.20 in 
which the axial velocity is represented by isovelocity contours. We can observe 
that at <^=22.5° the contours have shifted towards the outer bend becoming 
elliptic in shape. Further downstream we can see that the mmimum of axial 
velocity has shifted towards the outer wall, while a region of relatively low axial 
velocity has appeared near the inner wall. At <^=40.5° the region of low axial 
velocity has expanded from the inner wall along the plane of symmetry towards 
the tube. At the last station (f>=82.5° near the inner wall the axial velocity has 
increased. The results are shown in comparison with the experimental[15] in 
the fig.3.21. 

The secondary flow in the bend is shown in the fig.3.4.2 as velocity vec- 
tors in the cross plane. At (^=4.5° the vector plot shows the incipience of a 
vortex pair which gets intensified further downstream. At c!>=22.5° the vortex 
has shifted outwards. At <^=58.5° the higher secondary velocities appear at the 
inner half of the bend. At 4>=82.5° the vortex has developed a ‘tail' towards 
the upper wall (as seen by[15]) and the secondary flow has sharply reduced but 
remains complex in structure. The figures are substantially similar in detail 
to the experimental data of [15]. 
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Figure 3.21: Comparison of 
and experimental[15]. The i 










Chapter 4 


Conclusion 


The general form finite volume Navier Stokes solver has been thoroughly tested 
for laminar flows in 2D and 3D regular geometries. The results for 3D flows 
in complex geometries are close to experimental results of established repute. 
Thus this thesis validates the accuracy and robustness of the Navier Stokes 
solver. Thus the solver has now been fully validated and may be used to study 
laminar and turbulent flows in complex flows. 
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